/**********************************************************************
Compile : gcc -O3 -o generate6 check.c generate6.c
Usage   : ./generate6 mode mode1 mode2 pnflag N L M n seed       
mode    : 0 common; 1 cycled; 2: triggled cycle
mode1   : 0 disallow self-activation; 1 allow self-activation
mode2   : Additional Requirements: 0 no requirements; 1 restriction I; 2 restriction II; 3 restriction III
pnflag: 1 print the candidate processes
N       : number of the nodes                
L       : number of the steps        
M       : the number of "1" in the process  
n       : the number of the attempts         
seed    : random seed
***********************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>

long int counter2 = 0;

void print_help ()
{
	printf ("Usage         : ./generate6 mode mode1 mode2 pnflag N L M n seed\n");
	printf ("                (./generate6 0 0 2 1 11 12 41 100000 $RANDOM)\n");
	printf ("mode          : 0 normal ; 1 cycled ; 2 triggled cycled\n");
	printf ("mode1         : 0 disallow self-activation ; 1 allow self_activation\n");
	printf ("mode2         : 0 no Restrictions ; 1 Restriction I ; 2 Restriction II; 3 Restriction I and II\n");
	printf ("pnflag        : 1 print the candidate processes\n");
	printf ("N             : number of the nodes\n");
	printf ("L             : number of the steps\n");
	printf ("M             : the number of \"1\" in the process\n");
	printf ("n             : the number of the attempts\n");
	printf ("seed          : random seed\n\n");
	printf ("Note          : This program is special optimized for sampling under restriction II\n");
}

/* if one state can be a middle state of two different position, return 1 */
int check_pathway1 (int N, int i, int pathway[i + 1][N])
{
	int j, k, flag;

	flag = 0;

	for (j = 0; j < i - 2; j++)
	{
		flag = 1;
		for (k = 0; k < N; k++)
		{
			if (pathway[j][k] == pathway[j + 1][k] && pathway[i][k] == pathway[i - 1][k] && pathway[j][k] != pathway[i][k])
			{
				flag = 0;
				break;
			}
		}
		if (flag == 1)
			break;
	}

	return flag;
}

/* check all-0 or all-1 state */
int check_pathway2 (int N, int L, int pathway[L][N])
{
	int i, j, flag;

	for (i = 0; i < N; i++)
	{
		flag = 0;
		for (j = 0; j < L; j++)
		{
			flag += pathway[j][i];
		}
		if (flag == 0 || flag == L)
			break;
	}

	if (flag == 0 || flag == L)
		return 1;
	else
		return 0;
}

void generate_pathway2 (int N, int L, int M, int pathway[L][N], int mode)
{
	int i, j, n1, n2, flag;
	float ratio;

	do
	{
		counter2 ++;
		n1 = N * (L - 1);
		n2 = M;
		for (i = 0; i < L - 1; i++)
		{
			for (j = 0; j < N; j++)
			{
				ratio = (float) n2 / (float) n1;
				if (drand48 () < ratio)
				{
					pathway[i][j] = 1;
					n2 --;
				}
				else
					pathway[i][j] = 0;
				n1 --;
			}

			flag = check_pathway1 (N, i, pathway);
			if (flag == 1)
				break;
		}

		if (flag == 0)
		{
			flag = check_pathway2 (N, L - 1, pathway);
			if (flag == 1)
				counter2 --;
		}
	}
	while (flag);

	if (mode == 0)
	{
		for (i = 0; i < N; i++)
		{
			pathway[L - 1][i] = pathway[L - 2][i];
		}
	}
	else if (mode == 1)
	{
		for (i = 0; i < N; i++)
		{
			pathway[L - 1][i] = pathway[0][i];
		}
	}
}

void generate_pathway (int N, int L, int M, int pathway[L][N], int mode)
{
	if (mode != 2)
		generate_pathway2 (N, L, M, pathway, mode);
	else
	{
		int i, j;
		int pathway2[L - 1][N - 1];

		generate_pathway2 (N - 1, L - 1, M, pathway2, 0);

		pathway[0][0] = 1;
		for (i = 1; i < L; i++)
		{
			pathway[i][0] = 0;
		}

		for (i = 0; i < L - 1; i++)
		{
			for (j = 0; j < N - 1; j++)
			{
				pathway[i + 1][j + 1] = pathway2[i][j];
			}
		}

		for (i = 1; i < N; i++)
		{
			pathway[0][i] = pathway2[L - 2][i - 1];
		}
	}
}

int check_duplicate (int N, int L, int pathway[L][N])
{
	int i, j, flag;

	for (i = 0; i < L - 2; i++)
	{
		flag = 1;
		for (j = 0; j < N; j++)
		{
			if (pathway[i][j] != pathway[i + 1][j])
			{
				flag = 0;
				break;
			}
		}
		if (flag == 1) break;
	}

	return flag;
}

void combine_dist2 (int N1, int N2, int dist[N2 + 1], double Dist[N1 * N2 + 1], int i)
{
	double Dist2[N1 * N2 + 1];
	int j, k;

	if (i == 0)
	{
		for (j = 0; j <= N1 * N2; j++)
		{
			Dist[j] = 0.0;
		}
		for (j = 0; j <= N2; j++)
		{
			Dist[j] = dist[j];
		}
	}
	else
	{
		for (j = 0; j <= (i + 1) * N2; j++)
		{
			Dist2[j] = 0.0;
			if (j <= N2)
			{
				for (k = 0; k <= j; k++)
				{
					Dist2[j] += dist[k] * Dist[j - k];
				}
			}
			else
			{
				for (k = 0; k <= N2; k++)
				{
					Dist2[j] += dist[k] * Dist[j - k];
				}
			}
		}
		for (j = 0; j <= (i + 1) * N2; j++)
		{
			Dist[j] = Dist2[j];
		}
	}
}

void setup_count (int N, int L, double Count[N * (L - 1) + 1])
{
	int i, j, count[L];

	count[0] = count[1] = 1;

	for (i = 1; i < L - 1; i++)
	{
		count[i + 1] = 1;
		for (j = i; j > 0; j--)
		{
			count[j] = count[j] + count[j - 1];
		}
	}

	count[0] = count[L - 1] = 0;

	for (i = 0; i < N; i++)
	{
		combine_dist2 (N, L - 1, count, Count, i);
	}
}

int main (int argc, char *argv[])
{
	long int seed, n, counter;
	int N, L, M, i, j, k, p;
	int mode, mode1, mode2, pnflag;

	if (argc == 1)
	{
		print_help ();
		exit (0);
	}

	mode = atoi (argv[1]);
	mode1 = atoi (argv[2]);
	mode2 = atoi (argv[3]);
	pnflag = atoi (argv[4]);
	N = atoi (argv[5]);
	L = atoi (argv[6]);
	M = atoi (argv[7]);
	n = atoi (argv[8]);
	seed = atoi (argv[9]);
	srand48 (seed);
	printf ("Random Seed : %ld\tAttempts : %ld\n", seed, n);

	int pathway[L][N];

	counter = 0;

	for (i = 0; i < n; i ++)
	{
		generate_pathway (N, L, M, pathway, mode);

		if (check_duplicate (N, L, pathway) == 1)
			continue;

		p = network_inferred (N, L, pathway, mode, mode1, mode2);

		if (p == 1)
		{
			if (pnflag == 1)
			{
				printf ("c candidate process %d\n", M);
                                
				for (j = 0; j < L; j++)
				{
					for (k = 0; k < N; k++)
					{
						printf ("%d ", pathway[j][k]);
					}
					printf ("\n");
				}
			}
			
			counter ++;
		}
	}

	if (mode == 2)
	{
		N --;
		L --;
	}

	double Count[N * (L - 1) + 1];

	setup_count (N, L, Count);

	printf ("Total %ld candidate processes with activity %3d are obtained.\n", counter, M);
	printf ("Estimated feasible processes with activity %3d : %lg\n", M, Count[M] / (double) counter2 * (double) counter);
	printf ("Estimated possibility : %lg %lg\n", (double) counter / (double) counter2, (double) counter2);

	exit (0);
}
